Disease extinction in the presence of non-Gaussian noise 



OO 

o 
o 

(N 



Oh! 
I 

o 



C/3 

O 
• I— I 



> 
O 

O 

00 

O 



X 



Mark Dykman^, Ira B. Schwartz^, Alexandra S. Landsman^ 
^Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824 
^ US Naval Research Laboratory, Code 6792, Nonlinear Systems Dynamics Section, 
Plasma Physics Division, Washington, DC 20375 

We investigate stochastic extinction in an epidemic model and the impact of random vaccinations 
in large populations. We show that, in the absence of vaccinations, the effective entropic barrier for 
extinction displays scaling with the distance to the bifurcation point, with an unusual critical expo- 
nent. Even a comparatively weak Poisson-distributed vaccination leads to an exponential increase 
in the extinction rate, with the exponent that strongly depends on the vaccination parameters. 
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Practically all diseases of interest exhibit randomness. 
Childhood diseases [H, d, [1] , meningitis [1] , dengue fever 
and malaria 0] are but a few examples where incidence 
rates fluctuate. These fluctuations arise from fluctua- 
tions in population, epidemic parameters, and intrinsi- 
cally random contacts within the population @, 01 ■ As 
diseases evolve in large populations, there is the possibil- 
ity of extinction and reintroduction of the disease [1, . 
Extinction occurs where the number of infectives become 
so small that there is insufficient transmission to keep the 
disease in its endemic state [l3,[ll|,[ll. To gain insight 
into disease fade-out, one can think of the dynamics as 
coming from a nonlinear physical system far from ther- 
mal equilibrium, with extinction resulting from a large 
infrequent fluctuation. 

A well-known model of population dynamics is the 
so-called SIS model where only susceptibles (S) and in- 
fectives (I) are present P, [l^l- Here, in the absence 
of fluctuations the disease spread is characterized by 
the reproductive rate of infection, Rq, which is defined 
so that an endemic state exists along with the disease 
free equilibrium for Rq >1. The disease free state 
is unstable. However, in the presence of fluctuations, 
this state may be reached, albeit for a limited time 

H, [3, M, M, [13, H E, Eoi 

A major characteristic of fluctuation-induced extinc- 
tion is the extinction rate, or the reciprocal mean first 
time the number of infectives approaches zero. It has 
been studied in the continuous limit using the Langevin 
approach with fluctuations induced by noise. A discrete 
birth-death SIS model was investigated recently and a 
comparison with the continuous model was performed 
by Doering et al. [1^. The analysis referred to the one- 
variable model with detailed balance, which allows one 
to obtain an explicit solution. However, it does not re- 
veal some generic features of the SIS model, including 
the scaling behavior of the extinction rate. 

The goal of this letter is two-fold: (i) to show that the 
extinction rate scales with the control parameter Rq and 
to find the corresponding scaling exponent, and (ii) to 
study the effect of vaccination on extinction rate. We 
will be interested in the important case where the vac- 



cine schedule is a Poisson process. As wc show, even 
comparatively weak vaccination can increase the extinc- 
tion rate exponentially strongly. This is a dynamical ef- 
fect which happens, because an appropriate sequence of 
vaccine pulses can "resonate" with the dynamics of the 
system followed during extinction. 

We consider a model where susceptibles (S) are born 
at rate ^, both susceptibles and infectives (I) die at the 
same rate /.i, and infectives recover at rate >r and im- 
mediately become susceptible. If susceptibles contact 
infectives, they may become infected at rate (3. Time- 
dependent vaccination reduces the number of suscepti- 
bles at rate S,{t) Q. This rate will be assumed small, on 
average. The events of birth, death, and contact happen 
at random. They are transitions between the states of 
the system with different S and /. Therefore the quan- 
tity of interest is the probability p{S, I) to have given S 
and /. It is given by the master equation 



P(X) 



[W{X-r;r)p{X-v)-W{X;r)p{X)] (1) 



Here, we introduced vector X = {Xi,X2) with compo- 
nents Xi = S, X2 — I and vector r = (ri, with com- 
ponents ri and r2 showing, respectively, the increments 
in S and / in a single transition. The transition rates 
W{X, r) are 

W{X; (1, 0)) = N[fi - W{X; (-1,0)) = fiX^, 

W{X; (0, -1)) = fiX2, W{X; (1, -1)) = KX2, 
W{X;i-l,l)) = PX1X2/N, (2) 

where N is the scaling factor which we set equal to the 
average population, 1. 

For large S,I oc N fluctuations of S, I are small on 
average. If these fluctuations are disregarded, one arrives 
at the deterministic equations for the mean values of S, I 

Xi = N[^i - e(t)] - ^iXi + KX2 - /3X1X2/N, 
X2 = -{p + x)X2 + fiXiX2lN. (3) 

These are standard equations of the SIS model. In the 
absence of vaccination, ^(t) = 0, for Rq = /3/(/i + x) > 1 
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they have a stable endemic solution = Nxa with 
xiA = Rq^^ ^2A = 1 — Rq^- In addition, Eqs. ([3]) have an 
unstable stationary state (saddle point) X5 = A^x^ with 
xis = l,a^25 = 0. This state corresponds to extinction 
of infectives. 

For ^ 1 and for small vaccination rate the distribu- 
tion /3(X) has a peak at the stable state Xyi with width 
cx iV^/^. The probability of having a small number of 
infected, X2 -C X2A, is determined by the far tail of this 
peak. The tail can be obtained by seeking the solution 
of Eq. ([T]) in the eikonal form. 



p(X) exp[-A^.s(x)], X 
p(X + r) «p(X)exp(-pr), 



X/N, 
P = d^s. (4) 



For time- independent parameters W this formulation was 
used in a number of papers [11, ^ [U, [H, [H, [H [H] . 
However, in the present case it has special features, as 
shown below. 

To leading order in N~^, the equation for s has a 
form of the Hamilton-Jacobi equation s = — i7(x, 9xS; t), 
where s is the effective action, and the effective Hamil- 
tonian is 



i/(x,p;t)=^«;(x;r) (eP"" - 1) , 



(5) 



with w(x; r) = N~^W(X.;r) being the transition rates 
per person. Action s(x) can be found from classical tra- 
jectories of the auxiliary system with Hamiltonian H that 
satisfy equations 



i = dpH, 



P = -d^H. 



(6) 



We start with the case where there is no vaccine, 
^(t) = 0. In this case the transition rates w = and 
the Hamiltonian H — ij'"^ are independent of time. Of 
interest to us is the stationary distribution. The func- 
tion s = s^^^ is independent of time. It has the form 



= (0) 



pxdt, i7(")(x,p) = 0. (7) 



Here, the integral is calculated for a Hamiltonian trajec- 
tory (x(<), p{t)) that starts at t ^ —00 at x ^ x^i, p 
and arrives at time tf at a state JCf. This trajectory de- 
scribes the most probable sequence of elementary events 
X X-l-r bringing the system to X/ = A''x/. It provides 
the absolute minimum to s'-''-'(x/) [s^°^(x/) is indepen- 
dent oftf]. The quantity Ns'^^^x) is the entropic barrier 
for reaching X = A^x; it gives also the exponent in the 
expression for the mean first passage time for reaching X 
from the vicinity of the attractor X^ [2^ . 

The extinction rate is determined by s^^-* for X2 0. It 
is intuitively clear and can be shown from Eq. ^ that the 
minimum of s'*''(x) over xi for 2:2 —> is reached at the 



saddle point X5 of the fluctuation- free motion Thus 
the entropic barrier for extinction is A^Soxt — A^.s'°-'(x5). 

The Hamiltonian trajectory x^xt (i),Pcxt(0 that gives 
s^°-'(x5) is the optimal extinction trajectory. One can 
show that it approaches X5 as < — *■ 00. This is similar to 
other problems of an optimal trajectory leading from a 
deterministic stable state to a saddle point [2J,|27[. How- 
ever, in distinction from the more common situation, for 
t 00 the momentum pcxt does not go to zero. Instead 
Pcxt(i) P5, with ps = (0, — log Rq) . This is in spite of 
the fact that, along with {xs,ps), the Hamiltonian ij'"' 
has a "standard" fixed point (xg, p = 0). 

Indeed, the optimal extinction trajectory should lie on 
the stable manifold of the fixed point with x = X5 . The 
stable manifold of (x5,p = 0) lies on the plane X2 = 
Pi = 0. The point (x4,p = 0), and thus the optimal 
extinction trajectory do not lie on this plane. Therefore 
this trajectory may not go to (xg, p = 0). In contrast, it 
may go to the fixed point (xg, ps), whose stable manifold 
is not confined to a plane in the (x, p) space. 

The situation where an auxiliary Hamiltonian system 
has two fixed points with the same x^ was first noticed for 
a system described by the Fokker-Planck equation with 
a singular at X5 diffusion matrix and the "right" 
point was chosen based on numerical simulations. This 
situation was also found for a system described by a 
one- variable master equation, where the Hamiltonian dy- 
namics is integrable |18j : it occurs also in a two- variable 
susceptible-infected-recovered (SIR) inodel concurrently 
studied by Kamenev and Meerson [28| . 

Equations @ allow finding the extinction rate for any 
values of the parameters of the system. An explicit ana- 
lytical solution in the absence of vaccination can be ob- 
tained near the saddle-node bifurcation point Rq = 1 
where states xa and X5 merge. For 77 = /3(i?o — ^)/Ro ^ 
1 we have X2A = v/P 1- The relaxation time of X2 
is 1]^^ . It is much longer than the relaxation time of xi, 
which is /i^^, i.e., X2 is a soft mode, and xi follows X2 
adiabatically. 

In the adiabatic approximation we have in Hamiltonian 
equations ^ xi = 1 — X2,Pi = /3x2P2/lJ', while X2 ^ 
1, IP2I ^ 1- The equations for slow variables X2,P2 have 
a Hamiltonian form 



±2 = dH^''/dp2, P2 = -dH^''/dx2 



ad , 



(8) 



with Hamiltonian iJ*"^ = i]X2P2 ~ Px2P2{x2 ~ P2)- The 
Hamiltonian trajectory is 

P2{t)=X2{t)^'j, X2(t)=X2A(l + e''(*-*°))"\ (9) 

From Eqs. 0, ® 



„(0) _ J2 icyo2 



7772/3^ = (i?o - l)72i?o. 



(10) 



The entropic barrier for extinction Ns\,J^ pI7|) scales with 
the distance to the bifurcation point ?/ oc i?o — 1 as ry^. 
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This is in contrast to the standard scahng of the acti- 
vation energy of escape near a saddle-node bifurcation 
point, where the critical exponent is 3/2 [23|, as has 
been seen in various dynamical systems close and far from 
thermal equilibrium. Such unusual scaling is related to 
Ps being nonzero. It emerges also in the SIR model [2^. 

Generally, the scaled barrier Sg°|. depends on two pa- 
rameters, Rq and /i' = -I- >c). In Fig. [1] asymptotic 
expression (jlOp is compared with the results obtained 
from Eqs. ([6]), ([7]) for several /i'. There is a reasonably 
good agreement even far from Rq = 1. As Rq increases 
the dependence of Sg^i on ^' becomes more pronounced. 



Equations (g]), ^ give 

(p(X)) = A(X)pW(X), AiXf) = VSNXfit)], (13) 

where is the distribution for £_{t) = and 'P^[K{t)] = 
(exp [i J K{t)£,{t)dt\ ) is the characteristic functional of 
£_{t). Because iV ^ 1, an already weak noise ^{t) can 
significantly change the distribution (p(X)), the factor A 
can be exponentially large. 

To extend the analysis to the problem of extinction, we 
choose the final point to be X5. The corresponding loga- 
rithmic susceptibility Xoxt(i) is determined by the unper- 
turbed optimal extinction path (Xoxt(i), Poxt(O)- Since 
X5 is approached along the optimal path asymptotically 
as t — > cx), the action Scxt for reaching this point is given 
by Eq. (|12[) in which integration over time goes from —00 
to 00. The noise-induced change of the extinction rate is 
then determined by the factor 

Acxt = V^[iNXeAt)], Xoxt(t) = /l(xoxt(i), Pcxt(t))(14) 

We call Acxt the noise-induced extinction factor. 

The effect of noise on extinction can be illustrated us- 
ing an important model where the noise is a Poisson pro- 
cess: a sequence of pulses g{t — tj) occurring at random 
FIG. 1: The scaled barrier for extinction 5^°]. vs the reproduc- times tj. We assume that the pulses have a small ampli- 



(0) 

%x\ 0.2 
0.1 
0.0 



0.3|- x-^'=0.8 
□ -^'=0.6 



A - ^'=0.2 
y 

A 



2 



^0 



tive rate of infection 7?o. The dashed line shows asymptotic 
expression (fTO)) . The data points for different /x' = ^,/{iJ, + >c) 
are obtained from a numerical solution of Eqs. (|6]), ((Tjl. 

We now discuss the effect of vaccination. We are in- 
terested in the distribution (p(X)) averaged over realiza- 
tions of noise ^(i). If ^(t) is a stationary noise, (/3(X)) is 
stationary, too. The mean first time of reaching x from 
the vicinity of the stable state is assumed to largely ex- 
ceed the correlation time tcorr of ^ {t) . 

The full Hamiltonian of the system that determines s 
can be written as H = H'^^^ + with 

7j(i'(x,P,i) = -CWMx,p), = cxp(pi) - 1. (11) 

The term H^^^ is small for weak noise ^(t). Because 
s(xy,tj^) provides a minimum to the integral over time 
of p X — iJ , to first order in ^(t) we have [29[ 



tude and duration small compared to the relaxation time 
of the system and the reciprocal average pulse frequency 
;/~^ (the noise realizations of interest should satisfy the 
restriction ^{t) < fi). With account taken of the explicit 
form of the characteristic functional 3l| , 



^oxt — ^av^fli ^av — 6Xp 

/>oc 

Afi = exp 



1/ / dt ti{t) 

— 00 

dt (e-'"^*) - 1 + K{t)^ 



, (15) 



s(x/,t/)«sW(x;)+ / dtmXfit) 



(12) 



Here, Xf{t) = /i(x(t|x/, t/), p(t|x/, t/)) ; the function 
x(i|xy, i/), p(t|x/t/) describes the Hamiltonian trajec- 
tory © to X/ calculated for H = H^"K 

From Eq. (fT2)) . the logarithm of the distribution 
p(X.f,tf) is linear in the force ^{t). The proportional- 
ity coefficient is cx X/(Oj ^'^d therefore we call x/(t) the 
logarithmic susceptibility, as for systems where fluctua- 
tions are induced by Gaussian noise (soj ; it is convenient 
to set Xf{t) = for t > tf. 



with K{t) = NJ dt'xext{t)git - t') « Nxext{t)g, where 
g = J g{t)dt. The factor describes the effect of the 
average noise {^{t)), whereas the term Aa describes the 
effect of the fluctuating part of ^(i) with zero mean. 

The general expression (|15p is simplified in the limit- 
ing cases. For weak noise, where \K{t)\ <C 1, we have 
Afi — exp [1/ J d< K^(t)/2] , implying logAfl gc g^. In 
the opposite limit, max[— K(t)] ^ 1, we have Afi = 
exp [z/(27r/Km)^/^ exp(-Km)] , where = -K(t„i) is 

the maximum of —nit) and km = k{tm)- In this case 
log Afi is exponential in the pulse intensity g. Note that 
the results apply for logAext ^ Ns''°\ 

An explicit expression for A^xt can be obtained near 
the bifurcation point, 77 <C 1. From Eq. here Xext = 
±2//^. This gives A^^ = exp [1/775 A^//3/x]. The exponent in 
linearly scales with the distance to the bifurcation 
point rj and the pulse intensity g. 

For pulse duration small compared to the fluc- 
tuation part of the extinction factor Afi is determined 
by the parameter a = grf'N/^p. If cr ^ 1, then 
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Afi w exp[i/(T^/1277]. In the opposite limit, cr 1, we 
have Afi « exp[(4i//77)(7r/cr)^/^ exp(cr/4)]. Here the de- 
pendence of Afi on the vaccination pulse intensity g is 
double exponential. 

The logarithm of An as a function of a is shown in 
Fig. [21 The parabolic small-tr asymptotics works reason- 
ably well, numerically, for cr < 3, whereas the exponential 
large-(T asymptotics is approached for a > 15. 




FIG. 2: The scaled fluctuation-induced change of the barrier 
for extinction, Ah = (rj/i/) log An, vs the scaled intensity of 
noise pulses a — grfN/^f) (solid line). The results refer to 
r] <til. The dashed lines show the asymptotic behavior of Aa 
for small and large a. 

The exponentially strong effect of random vaccine on 
disease extinction results from a dynamical cooperation 
between an outburst of noise of an appropriate tempo- 
ral shape and the evolution of the system along the op- 
timal path leading to extinction. One can think of the 
vaccine- induced change of the extinction barrier as a gen- 
eralized work done by the corresponding noise realization 
along the optimal path. In this picture, a zero-mean noise 
should be expected to reduce the extinction barrier, and 
then the respective part of the noise-induced extinction 
factor Aff should exceed 1. This is indeed the case for 
Poisson noise, as seen from Eq. 

In summary, we have considered fluctuations in the full 
two- variable SIS model and found the rate of extinction 
of disease with and without vaccination. We have devel- 
oped a general approach in which the problem is reduced 
to the analysis of dynamics of an auxiliary Hamiltonian 
system, with nontrivial boundary conditions. We show 
that the entropic barrier for extinction displays scaling 
dependence on the reproductive rate of infection i?o for 
small i?o — Ij with an unusual exponent. Even compara- 
tively weak vaccination can exponentially strongly affect 
the extinction rate. A general expression that describes 
this effect for random vaccine in terms of its characteristic 
functional has been obtained. For a Poisson distributed 
vaccine, the change of the exponent of the extinction rate 
may itself depend on the vaccine strength exponentially. 
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